How to Integrate a Derivative

In the previous post, I had tried to use derivatives to improve the estimate for a single sample, but that turned out not to work very well. That was an approximation. To use the derivative exactly, the definition is f(x) - f(0) = int_0^x f`(y) dy.

Does this have a practical application? After all, if we do this on the example in the previous post, we go from a single integral to a double integral - it looks like we're making the problem harder! Ah, but this is Monte Carlo. The convergence rate is independent of dimensionality. So adding more integrals doesn't necessarily make the problem harder.

The final error does, however, depend on the intrinsic variance of the integral. In this example, this variance gets much worse, so making this transformation is not a good idea. In general, derivatives will vary more than the original function, so this technique is not likely useful in general. Unless the transformed problem has some other desirable features, like better importance sampling or something.

Finally, I should note that something similar to this is used in thermodynamic integration to get free energies.

Sometimes it doesn't work

The other day I had an idea about trying to improve MC samples. The Monte Carlo method for evaluating an integral is the average of a set of samples. Each sample can be considered a very crude approximation to the integral - each sample approximates the integrand as a constant. What if we could make each sample a better approximation, say by using a the derivative information at the sample point?

Another consideration is that the intrinsic variance of the integral (int [f(x)-ave(f)]^2) controls the accuracy of the MC approximation. If we could smooth the integrand, and lower the variance, the MC integral would be better. Perhaps expanding the integral in a series around the sample point and integrating would help smooth it?

Then reality struck. I tried it on v(x)*exp(-v(x)/T), with v(x) = 1/x^12. I made the expansion and integrated in an interval from [x-h,x+h], and graphed the result - it was clearly *not* smoother than the original function. It was also clear that the expansion was an approximation, and only correct in the limit of small h.

Upon further reflection, I realized a few things. One is that each MC sample may be a crude approximation to the integral, but it has the important property of being unbiased (ie, all the errors will average to zero). Any improvement in the sample value needs to retain this property.

Secondly, smoothing would be most valuable when the integrand is rapidly changing. And this is just the condition where the Taylor expansion is likely to be the least accurate. I suppose that any local smoothing of the integrand will suffer from this sort of problem.

Sometimes, it doesn't work. Ah, but that's what research is all about.

Visualizing Wavefunctions

I've started learning VTK, Python, and Tk in order to look into visualizing wavefunctions. So far, I've made an application that plots a function of x and y. The user can change the plotted function in an edit box. Get the 2D plotter here.

I've also got a prototype application that plots a function of 3 variables, and lets the user choose a slice (a plane) through the displayed volume, and the application plots the 2D slice in a view that looks very much like the 2D plotter. But it's not quite ready for public consumption yet.

Welcome to the Metropolis Disco - How do you move your particles?

Out on the dance floor tonight we have the spotlight on two papers on improving the efficiency of Monte Carlo moves. But first, a brief summary of the results of Peskun's theorem (Both papers mention it.). The Metropolis algorithm takes a trial move and accepts or rejects the move with some probability. If the move is rejected, that location is used in the average again. Rejected moves will increase the autocorrelation, and so increase the variance of the resulting average. Peskun's theorem validates this intuition, under certain conditions. So the result is that, just like dancing, moving is better than standing still.

The first paper is Delayed Rejection Variational Monte Carlo by Bressanini, Morosi, and Tarasco. The problem is that for best efficiency, moves near an atomic nucleus should have small step, and moves far away from a nucleus should have a large step. If a move is rejected, it is likely that the step size was too large. One could imagine a dynamic or adaptive step size method that lowers the step size when there are lots of rejections and increases it when there are lots of acceptances (*). This paper does something slightly different - when a move is rejected, it is retried with a smaller step, but the original rejected move is not made part of the average. The acceptance probability is adjusted to take this into account. So there is an increase in the amount of time it takes for a move, but that is made up for by a larger decrease in the autocorrelation time.

The second paper is Optimal Monte Carlo Updating by Pollet, Rombouts, Van Houcke, and Heyde. They work with discrete systems, and introduce a a procedure called "Metropolizing" - which is to take a move rule (eg, heat bath) and trasform it to a new, more efficient transition matrix. This process can be applied repeated until there is only one non-zero diagonal element, resulting in a very efficient transition matrix (For single site updating rules).



(*) I'm not entirely sure a simple dynamic step size method would satisfy detailed balance. I think it would, but that should be checked.

More Graphic Card Computations

There's a site, GPGPU (General Purpose computation on Graphics Processing Units), dedicated to interesting things one can do with a graphics card. (Well, there's always graphics, but using things for their intended purpose can get boring :)

Last fall, I wondered about the possiblity of doing Monte Carlo on the graphics card.

On the GPGPU site there's a link to a paper about doing Ising model and percolation simulations on the graphics card (Benchmarking and Implementation of Probability-Based Simulations on Programmable Graphics Cards). The work is primarily for proof of concept and benchmarking, with an eye towards lattice QCD simulations eventually. Hmm. Play a killer game of Quake *and* uncover fundamental secrets of the universe - what a combination!

I can imagine students requesting the latest "Ising co-processors" so they can "make progress on their thesis" when Doom 3 comes out. (Err, maybe they need to leave off that last part when making the request).

PIMC live

Check out the Path Integral simulation applet on John Shumway's page. You can change the temperature and see the paths localize at high temperature and spread out at low temperature.

Pseudorandom number generators and the stock market

As correlations are discovered in psuedorandom number generators used in simulations, the correlations are removed by creating better generators. So as time marches on, the state of the art PRNG's should be getting more and more random.

As correlations are discovered in the stock market, they are removed as the resevoir of correlation is drained for its profit. So the stock market should be getting more and more random as time goes on (especially as faster computers and more sophisticated statistical techniques are used)